I, \[**insert your name**\], confirm that the work presented in this assessment is my own. Where information has been derived from other sources, I confirm that this has been indicated in the work.
date: 16 December, 2021
# load the package
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.3-0
## Loading required package: spatstat.core
## Loading required package: nlme
## Loading required package: rpart
## spatstat.core 2.3-1
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-0
##
## spatstat 2.2-0 (nickname: 'That's not important right now')
## For an introduction to spatstat, type 'beginner'
library(here)
## here() starts at /Users/mwanghoffin/Documents/GitHub/CASA0005/Practice Exam 1
library(sp)
library(rgeos)
## rgeos version: 0.5-8, (SVN revision 679)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.4-5
## Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
library(GISTools)
## Loading required package: RColorBrewer
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:spatstat.geom':
##
## area
library(tmap)
## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sf
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.3, PROJ 7.2.1
library(geojson)
##
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
##
## polygon
library(geojsonio)
## Registered S3 method overwritten by 'geojsonsf':
## method from
## print.geojson geojson
##
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
##
## pretty
library(tmaptools)
# first, get the Boundaries of Community Districts in New York City
NYDistricts <- st_read(here::here("Community Districts", "geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b.shp")) %>%
st_transform(., 32618)
## Reading layer `geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b' from data source `/Users/mwanghoffin/Documents/GitHub/CASA0005/Practice Exam 1/Community Districts/geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 71 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS: WGS84(DD)
NYDistricts
## Simple feature collection with 71 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 563069.7 ymin: 4483098 xmax: 609762.3 ymax: 4529952
## Projected CRS: WGS 84 / UTM zone 18N
## First 10 features:
## boro_cd shape_area shape_leng geometry
## 1 206 42664312 35875.71 MULTIPOLYGON (((595104.1 45...
## 2 404 65739662 37018.37 MULTIPOLYGON (((597308.6 45...
## 3 203 44796871 33500.07 MULTIPOLYGON (((594364.8 45...
## 4 304 56662613 37007.81 MULTIPOLYGON (((593253.9 45...
## 5 205 38316975 29443.05 MULTIPOLYGON (((593432.1 45...
## 6 312 99525501 52245.83 MULTIPOLYGON (((586966.1 45...
## 7 483 191997989 106865.65 MULTIPOLYGON (((605960.8 44...
## 8 411 260363028 103434.76 MULTIPOLYGON (((605933.1 45...
## 9 207 53311689 44812.15 MULTIPOLYGON (((594782.6 45...
## 10 208 92071724 47817.43 MULTIPOLYGON (((592919.8 45...
qtm(NYDistricts)
# Now lets get the location of the evictions
library(tidyverse)
## Registered S3 method overwritten by 'cli':
## method from
## print.boxx spatstat.geom
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.5 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::collapse() masks nlme::collapse()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select()
library(janitor)
##
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
Evictions <- read_csv(here::here("Evictions.csv"),
na = c("", "NA", "n/a"),
locale = locale(encoding = 'Latin1'),
col_names = TRUE) %>%
na.omit(Evictions)
## Rows: 66559 Columns: 20
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (12): Court Index Number, Docket Number, Eviction Address, Eviction Apar...
## dbl (8): Eviction Postcode, Latitude, Longitude, Community Board, Council D...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#NewEvictions <- na.omit(Evictions)
NewEvictions <- Evictions %>%
st_as_sf(., coords = c("Longitude","Latitude"),
crs = 4326) %>%
st_transform(., 32618) %>%
clean_names()
summary(NewEvictions)
## court_index_number docket_number eviction_address
## Length:52917 Length:52917 Length:52917
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## eviction_apartment_number executed_date marshal_first_name
## Length:52917 Length:52917 Length:52917
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## marshal_last_name residential_commercial borough eviction_postcode
## Length:52917 Length:52917 Length:52917 Min. :10000
## Class :character Class :character Class :character 1st Qu.:10453
## Mode :character Mode :character Mode :character Median :10468
## Mean :10739
## 3rd Qu.:11221
## Max. :11698
## ejectment eviction_legal_possession community_board council_district
## Length:52917 Length:52917 Min. : 1.000 Min. : 1.00
## Class :character Class :character 1st Qu.: 4.000 1st Qu.:12.00
## Mode :character Mode :character Median : 7.000 Median :17.00
## Mean : 7.899 Mean :22.91
## 3rd Qu.:12.000 3rd Qu.:36.00
## Max. :18.000 Max. :51.00
## census_tract bin bbl nta
## Min. : 1 Min. :1000000 Min. :1.000e+09 Length:52917
## 1st Qu.: 195 1st Qu.:2008102 1st Qu.:2.027e+09 Class :character
## Median : 371 Median :2100990 Median :2.046e+09 Mode :character
## Mean : 8546 Mean :2578910 Mean :2.508e+09
## 3rd Qu.: 988 3rd Qu.:3216952 3rd Qu.:3.051e+09
## Max. :157901 Max. :5171959 Max. :5.080e+09
## geometry
## POINT :52917
## epsg:32618 : 0
## +proj=utm ...: 0
##
##
##
#plot the evictions in the city
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(NYDistricts) +
tm_polygons(col = NA, alpha = 0.5) +
tm_shape(NewEvictions) +
tm_dots(col = "blue")
This part of R might not useful.
# join the New York District data and eviction data
#Joinfun <- function(data1, data2){
#output<- data1%>%
# st_join(NYDistricts,.)%>%
# add_count(boro_cd, name="eviction-borough")
# return(output)
#}
#NYEvictions <- Joinfun(NewEvictions, NYDistricts)
#NYEvictions
This two part cannot be loaded.
#remove duplicates
#library(tidyverse)
#library(sf)
#NYEvictions <- distinct(NYEvictions)
#NYEvictionsSub <- NYEvictions[NYDistricts,]
#check to see that they've been removed
#tmap_mode("plot")
#tm_shape(NYDistricts) +
# tm_polygons(col = NA, alpha = 0.5) +
#tm_shape(NYEvictionsSub) +
# tm_dots(col = "blue")
Try only extract one district
#extract the district
CD502 <- NYDistricts %>%
filter(., boro_cd == 502)
#Check to see that the correct borough has been pulled out
tm_shape(CD502) +
tm_polygons(col = NA, alpha = 0.5)
summary(CD502)
## boro_cd shape_area shape_leng geometry
## Min. :502 Min. :591533049 Min. :142729 MULTIPOLYGON :1
## 1st Qu.:502 1st Qu.:591533049 1st Qu.:142729 epsg:32618 :0
## Median :502 Median :591533049 Median :142729 +proj=utm ...:0
## Mean :502 Mean :591533049 Mean :142729
## 3rd Qu.:502 3rd Qu.:591533049 3rd Qu.:142729
## Max. :502 Max. :591533049 Max. :142729
Clip the eviction data to the single borough
#clip the data to our single borough
NewEvictionsCD <- NewEvictions[CD502,]
#check that it's worked
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(CD502) +
tm_polygons(col = NA, alpha = 0.5) +
tm_shape(NewEvictionsCD) +
tm_dots(col = "blue")
The first thing we need to do is create an observation window for to carry out its analysis within — we’ll set this to the extent of the CD503 boundary
#now set a window as the borough boundary
window <- as.owin(CD502)
plot(window)
To use spatstat, we need point pattern analysis, we need to create a point pattern (ppp) object.
#create a sp object
NewEvictionsCD <- NewEvictionsCD %>%
as(., 'Spatial')
#create a ppp object
NewEvictionsCD.ppp <- ppp(x = NewEvictionsCD@coords[,1],
y = NewEvictionsCD@coords[,2],
window=window)
## Warning: data contain duplicated points
# just to check
NewEvictionsCD@coords[,1]
## [1] 573358.6 578156.8 568457.7 575043.7 576533.2 573475.7 571264.8 578739.8
## [9] 573445.4 579349.6 571217.3 571217.3 579368.4 570707.0 577399.5 574357.7
## [17] 579349.6 578940.5 577432.9 578638.3 576131.2 577828.3 578156.8 577999.9
## [25] 570766.8 571591.3 578483.1 579325.3 573570.1 577271.1 575865.3 575512.0
## [33] 574361.9 575330.0 568457.7 577218.0 579349.6 573191.5 571253.7 570630.6
## [41] 570895.6 576826.6 570950.7 578959.5 570865.5 576283.9 575507.1 577543.6
## [49] 579368.4 579205.0 578813.4 575512.0 578887.1 571622.2 579323.1 579090.7
## [57] 570613.7 579289.4 578156.8 571130.6 571396.9 576076.3 575889.1 576106.0
## [65] 576885.2 578461.4 570907.8 571468.0 570865.5 574357.7 578478.0 575053.3
## [73] 575743.3 570486.8 570613.7 574577.3 571026.6 575043.7 571372.9 576077.3
## [81] 573156.1 578156.8 576062.9 571889.8 578813.4 575878.2 574364.0 578941.7
## [89] 574433.5 571396.9 575507.1 576471.6 578910.9 576543.9 572858.7 578481.6
## [97] 571117.6 568559.0 575865.3 571832.7 571073.3 575865.3 576283.9 574364.0
## [105] 576077.4 568728.8 575707.5 577857.6 578057.7 574353.5 570970.1 577949.3
## [113] 575507.1 571199.3 576077.3 575512.0
# check the new ppp object
NewEvictionsCD.ppp %>%
plot(.,pch=16,cex=0.5,
main="Evictions in CD502")
The size and shape of the Kernel affects the density pattern produced
NewEvictionsCD.ppp %>%
density(., sigma = 500) %>%
plot()
K <- NewEvictionsCD.ppp %>%
Kest(., correction="border") %>%
plot()
Kval <- as.data.frame(Kest(NewEvictionsCD.ppp, correction = "border"))
Kval
## r theo border
## 1 0.000000 0.000000e+00 244878.6
## 2 4.026481 5.093323e+01 261203.8
## 3 8.052962 2.037329e+02 285691.7
## 4 12.079443 4.583991e+02 285691.7
## 5 16.105925 8.149317e+02 302016.9
## 6 20.132406 1.273331e+03 302016.9
## 7 24.158887 1.833596e+03 318342.2
## 8 28.185368 2.495728e+03 383643.1
## 9 32.211849 3.259727e+03 391805.7
## 10 36.238330 4.125592e+03 391805.7
## 11 40.264812 5.093323e+03 398679.5
## 12 44.291293 6.162921e+03 406985.4
## 13 48.317774 7.334385e+03 435725.0
## 14 52.344255 8.607716e+03 452483.6
## 15 56.370736 9.982913e+03 452483.6
## 16 60.397217 1.145998e+04 413177.0
## 17 64.423698 1.303891e+04 429186.0
## 18 68.450180 1.471970e+04 446633.9
## 19 72.476661 1.650237e+04 468923.1
## 20 76.503142 1.838690e+04 523029.6
## 21 80.529623 2.037329e+04 523029.6
## 22 84.556104 2.246155e+04 556978.8
## 23 88.582585 2.465168e+04 556978.8
## 24 92.609067 2.694368e+04 556978.8
## 25 96.635548 2.933754e+04 556978.8
## 26 100.662029 3.183327e+04 556978.8
## 27 104.688510 3.443086e+04 594110.7
## 28 108.714991 3.713032e+04 686940.5
## 29 112.741472 3.993165e+04 696223.4
## 30 116.767953 4.283485e+04 705506.4
## 31 120.794435 4.583991e+04 724072.4
## 32 124.820916 4.894683e+04 724072.4
## 33 128.847397 5.215563e+04 724072.4
## 34 132.873878 5.546629e+04 721866.5
## 35 136.900359 5.887881e+04 700284.8
## 36 140.926840 6.239321e+04 700284.8
## 37 144.953321 6.600947e+04 707656.2
## 38 148.979803 6.972759e+04 717623.2
## 39 153.006284 7.354758e+04 717623.2
## 40 157.032765 7.746944e+04 755476.5
## 41 161.059246 8.149317e+04 765549.5
## 42 165.085727 8.561876e+04 775622.5
## 43 169.112208 8.984622e+04 795768.6
## 44 173.138690 9.417554e+04 807921.9
## 45 177.165171 9.860673e+04 828505.9
## 46 181.191652 1.031398e+05 837610.4
## 47 185.218133 1.077747e+05 848015.5
## 48 189.244614 1.125115e+05 879230.8
## 49 193.271095 1.173502e+05 910041.4
## 50 197.297576 1.222907e+05 931082.8
## 51 201.324058 1.273331e+05 941603.5
## 52 205.350539 1.324773e+05 973165.7
## 53 209.377020 1.377235e+05 984100.1
## 54 213.403501 1.430714e+05 1021336.3
## 55 217.429982 1.485213e+05 979514.4
## 56 221.456463 1.540730e+05 1033931.8
## 57 225.482945 1.597266e+05 1033931.8
## 58 229.509426 1.654821e+05 1056964.3
## 59 233.535907 1.713394e+05 1101004.5
## 60 237.562388 1.772986e+05 1145044.7
## 61 241.588869 1.833596e+05 1145044.7
## 62 245.615350 1.895225e+05 1189084.9
## 63 249.641831 1.957873e+05 1244135.1
## 64 253.668313 2.021540e+05 1255145.2
## 65 257.694794 2.086225e+05 1277165.2
## 66 261.721275 2.151929e+05 1297760.5
## 67 265.747756 2.218651e+05 1297760.5
## 68 269.774237 2.286393e+05 1297760.5
## 69 273.800718 2.355153e+05 1347888.6
## 70 277.827200 2.424931e+05 1347888.6
## 71 281.853681 2.495728e+05 1347888.6
## 72 285.880162 2.567544e+05 1380367.8
## 73 289.906643 2.640379e+05 1344312.9
## 74 293.933124 2.714232e+05 1356002.6
## 75 297.959605 2.789104e+05 1356002.6
## 76 301.986086 2.864994e+05 1373537.1
## 77 306.012568 2.941903e+05 1408606.2
## 78 310.039049 3.019831e+05 1420295.8
## 79 314.065530 3.098778e+05 1443675.2
## 80 318.092011 3.178743e+05 1478744.2
## 81 322.118492 3.259727e+05 1525502.9
## 82 326.144973 3.341729e+05 1550489.6
## 83 330.171455 3.424750e+05 1550489.6
## 84 334.197936 3.508790e+05 1562325.4
## 85 338.224417 3.593849e+05 1597832.8
## 86 342.250898 3.679926e+05 1633340.2
## 87 346.277379 3.767022e+05 1602385.0
## 88 350.303860 3.855136e+05 1602385.0
## 89 354.330341 3.944269e+05 1626663.6
## 90 358.356823 4.034421e+05 1650942.2
## 91 362.383304 4.125592e+05 1687360.0
## 92 366.409785 4.217781e+05 1687360.0
## 93 370.436266 4.310989e+05 1699499.3
## 94 374.462747 4.405215e+05 1699499.3
## 95 378.489228 4.500460e+05 1735917.1
## 96 382.515709 4.596724e+05 1735917.1
## 97 386.542191 4.694006e+05 1740016.1
## 98 390.568672 4.792308e+05 1752313.0
## 99 394.595153 4.891627e+05 1752313.0
## 100 398.621634 4.991966e+05 1776906.9
## 101 402.648115 5.093323e+05 1801500.8
## 102 406.674596 5.195699e+05 1887579.3
## 103 410.701078 5.299093e+05 1992103.3
## 104 414.727559 5.403506e+05 2004400.2
## 105 418.754040 5.508938e+05 2016697.1
## 106 422.780521 5.615389e+05 2005856.4
## 107 426.807002 5.722858e+05 2043232.6
## 108 430.833483 5.831345e+05 2049462.0
## 109 434.859964 5.940852e+05 2086838.2
## 110 438.886446 6.051377e+05 2111755.7
## 111 442.912927 6.162921e+05 2124214.4
## 112 446.939408 6.275483e+05 2149131.9
## 113 450.965889 6.389064e+05 2139912.4
## 114 454.992370 6.503664e+05 2168830.1
## 115 459.018851 6.619283e+05 2181625.6
## 116 463.045333 6.735920e+05 2220012.0
## 117 467.071814 6.853575e+05 2230967.0
## 118 471.098295 6.972250e+05 2250423.1
## 119 475.124776 7.091943e+05 2250423.1
## 120 479.151257 7.212655e+05 2276364.6
## 121 483.177738 7.334385e+05 2276364.6
## 122 487.204219 7.457134e+05 2347155.6
## 123 491.230701 7.580902e+05 2413836.1
## 124 495.257182 7.705688e+05 2453844.5
## 125 499.283663 7.831493e+05 2467180.6
## 126 503.310144 7.958317e+05 2441556.2
## 127 507.336625 8.086160e+05 2476941.0
## 128 511.363106 8.215021e+05 2518109.0
## 129 515.389588 8.344900e+05 2545554.4
## 130 519.416069 8.475799e+05 2586722.4
## 131 523.442550 8.607716e+05 2638686.9
## 132 527.469031 8.740652e+05 2663938.0
## 133 531.495512 8.874606e+05 2706334.9
## 134 535.521993 9.009579e+05 2734599.4
## 135 539.548474 9.145571e+05 2734599.4
## 136 543.574956 9.282581e+05 2734599.4
## 137 547.601437 9.420610e+05 2741665.6
## 138 551.627918 9.559658e+05 2826245.3
## 139 555.654399 9.699724e+05 2854938.1
## 140 559.680880 9.840809e+05 2912323.8
## 141 563.707361 9.982913e+05 2912323.8
## 142 567.733843 1.012604e+06 2941016.6
## 143 571.760324 1.027018e+06 2984055.9
## 144 575.786805 1.041534e+06 2984055.9
## 145 579.813286 1.056151e+06 2998402.3
## 146 583.839767 1.070871e+06 3012748.7
## 147 587.866248 1.085693e+06 3166076.1
## 148 591.892729 1.100616e+06 3173473.5
## 149 595.919211 1.115641e+06 3188268.3
## 150 599.945692 1.130769e+06 3210460.4
## 151 603.972173 1.145998e+06 3225255.1
## 152 607.998654 1.161329e+06 3225255.1
## 153 612.025135 1.176761e+06 3236411.8
## 154 616.051616 1.192296e+06 3251934.2
## 155 620.078097 1.207932e+06 3251934.2
## 156 624.104579 1.223671e+06 3267456.5
## 157 628.131060 1.239511e+06 3282978.9
## 158 632.157541 1.255453e+06 3298501.3
## 159 636.184022 1.271497e+06 3360590.7
## 160 640.210503 1.287643e+06 3376113.1
## 161 644.236984 1.303891e+06 3391635.4
## 162 648.263466 1.320240e+06 3391635.4
## 163 652.289947 1.336692e+06 3391635.4
## 164 656.316428 1.353245e+06 3407157.8
## 165 660.342909 1.369900e+06 3506606.1
## 166 664.369390 1.386657e+06 3538703.2
## 167 668.395871 1.403516e+06 3538703.2
## 168 672.422352 1.420477e+06 3570800.3
## 169 676.448834 1.437539e+06 3586848.8
## 170 680.475315 1.454704e+06 3618945.9
## 171 684.501796 1.471970e+06 3683140.1
## 172 688.528277 1.489339e+06 3683140.1
## 173 692.554758 1.506809e+06 3699188.6
## 174 696.581239 1.524381e+06 3731285.7
## 175 700.607721 1.542054e+06 3731285.7
## 176 704.634202 1.559830e+06 3753639.0
## 177 708.660683 1.577708e+06 3787455.6
## 178 712.687164 1.595687e+06 3804363.8
## 179 716.713645 1.613768e+06 3855088.7
## 180 720.740126 1.631952e+06 3863542.8
## 181 724.766607 1.650237e+06 3871997.0
## 182 728.793089 1.668624e+06 3905813.6
## 183 732.819570 1.687112e+06 3939630.1
## 184 736.846051 1.705703e+06 3964992.5
## 185 740.872532 1.724395e+06 3964992.5
## 186 744.899013 1.743190e+06 3976828.3
## 187 748.925494 1.762086e+06 3994044.0
## 188 752.951976 1.781084e+06 4080122.6
## 189 756.978457 1.800184e+06 4080122.6
## 190 761.004938 1.819386e+06 4080122.6
## 191 765.031419 1.838690e+06 4131769.7
## 192 769.057900 1.858095e+06 4209240.4
## 193 773.084381 1.877603e+06 4235063.9
## 194 777.110862 1.897212e+06 4235063.9
## 195 781.137344 1.916923e+06 4286711.1
## 196 785.163825 1.936736e+06 4286711.1
## 197 789.190306 1.956651e+06 4381397.5
## 198 793.216787 1.976668e+06 4424436.7
## 199 797.243268 1.996786e+06 4424436.7
## 200 801.269749 2.017007e+06 4510515.3
## 201 805.296231 2.037329e+06 4587985.9
## 202 809.322712 2.057753e+06 4639633.1
## 203 813.349193 2.078279e+06 4717103.7
## 204 817.375674 2.098907e+06 4734319.5
## 205 821.402155 2.119637e+06 4734319.5
## 206 825.428636 2.140469e+06 4734319.5
## 207 829.455117 2.161403e+06 4760621.2
## 208 833.481599 2.182438e+06 4760621.2
## 209 837.508080 2.203575e+06 4786923.0
## 210 841.534561 2.224814e+06 4821992.0
## 211 845.561042 2.246155e+06 4874595.6
## 212 849.587523 2.267598e+06 4909664.6
## 213 853.614004 2.289143e+06 4909664.6
## 214 857.640485 2.310790e+06 4975501.8
## 215 861.666967 2.332538e+06 4975501.8
## 216 865.693448 2.354389e+06 5029097.8
## 217 869.719929 2.376341e+06 5029097.8
## 218 873.746410 2.398395e+06 5046963.2
## 219 877.772891 2.420551e+06 5082693.9
## 220 881.799372 2.442809e+06 5134915.7
## 221 885.825854 2.465168e+06 5217034.4
## 222 889.852335 2.487630e+06 5263449.3
## 223 893.878816 2.510193e+06 5291298.2
## 224 897.905297 2.532859e+06 5328430.1
## 225 901.931778 2.555626e+06 5365562.0
## 226 905.958259 2.578495e+06 5384128.0
## 227 909.984740 2.601466e+06 5439825.9
## 228 914.011222 2.624538e+06 5523372.7
## 229 918.037703 2.647713e+06 5560504.6
## 230 922.064184 2.670989e+06 5634768.4
## 231 926.090665 2.694368e+06 5634768.4
## 232 930.117146 2.717848e+06 5681183.3
## 233 934.143627 2.741430e+06 5785338.4
## 234 938.170109 2.765114e+06 5823212.9
## 235 942.196590 2.788900e+06 5851618.8
## 236 946.223071 2.812788e+06 5889493.4
## 237 950.249552 2.836777e+06 5955773.9
## 238 954.276033 2.860869e+06 6096644.0
## 239 958.302514 2.885062e+06 6096644.0
## 240 962.328995 2.909357e+06 6096644.0
## 241 966.355477 2.933754e+06 6125629.7
## 242 970.381958 2.958253e+06 6135291.5
## 243 974.408439 2.982854e+06 6173939.0
## 244 978.434920 3.007556e+06 6318867.2
## 245 982.461401 3.032361e+06 6328529.1
## 246 986.487882 3.057267e+06 6328529.1
## 247 990.514364 3.082275e+06 6338190.9
## 248 994.540845 3.107385e+06 6357514.7
## 249 998.567326 3.132597e+06 6396162.2
## 250 1002.593807 3.157911e+06 6415486.0
## 251 1006.620288 3.183327e+06 6415486.0
## 252 1010.646769 3.208844e+06 6434809.7
## 253 1014.673250 3.234464e+06 6560414.1
## 254 1018.699732 3.260185e+06 6608723.5
## 255 1022.726213 3.286008e+06 6628047.2
## 256 1026.752694 3.311933e+06 6666694.7
## 257 1030.779175 3.337960e+06 6618184.1
## 258 1034.805656 3.364089e+06 6618184.1
## 259 1038.832137 3.390319e+06 6618184.1
## 260 1042.858619 3.416652e+06 6647773.6
## 261 1046.885100 3.443086e+06 6706952.6
## 262 1050.911581 3.469623e+06 6668339.3
## 263 1054.938062 3.496261e+06 6728777.4
## 264 1058.964543 3.523001e+06 6728777.4
## 265 1062.991024 3.549842e+06 6728777.4
## 266 1067.017505 3.576786e+06 6748923.5
## 267 1071.043987 3.603832e+06 6789215.6
## 268 1075.070468 3.630979e+06 6849653.7
## 269 1079.096949 3.658228e+06 6900018.8
## 270 1083.123430 3.685579e+06 7029435.2
## 271 1087.149911 3.713032e+06 7101479.2
## 272 1091.176392 3.740587e+06 7204399.2
## 273 1095.202873 3.768244e+06 7245567.2
## 274 1099.229355 3.796003e+06 7531871.9
## 275 1103.255836 3.823863e+06 7618951.3
## 276 1107.282317 3.851825e+06 7618951.3
## 277 1111.308798 3.879890e+06 7805853.5
## 278 1115.335279 3.908056e+06 7955966.1
## 279 1119.361760 3.936324e+06 7979060.4
## 280 1123.388242 3.964694e+06 8036796.0
## 281 1127.414723 3.993165e+06 8059890.2
## 282 1131.441204 4.021739e+06 8059890.2
## 283 1135.467685 4.050414e+06 8059890.2
## 284 1139.494166 4.079191e+06 8082984.4
## 285 1143.520647 4.108071e+06 8036203.8
## 286 1147.547128 4.137052e+06 8036203.8
## 287 1151.573610 4.166134e+06 8036203.8
## 288 1155.600091 4.195319e+06 8110636.8
## 289 1159.626572 4.224606e+06 8110636.8
## 290 1163.653053 4.253994e+06 8222765.4
## 291 1167.679534 4.283485e+06 8222765.4
## 292 1171.706015 4.313077e+06 8222765.4
## 293 1175.732497 4.342771e+06 8322435.3
## 294 1179.758978 4.372567e+06 8322435.3
## 295 1183.785459 4.402465e+06 8372270.2
## 296 1187.811940 4.432464e+06 8397187.7
## 297 1191.838421 4.462566e+06 8447022.6
## 298 1195.864902 4.492769e+06 8447022.6
## 299 1199.891383 4.523075e+06 8470593.2
## 300 1203.917865 4.553482e+06 8470593.2
## 301 1207.944346 4.583991e+06 8521775.0
## 302 1211.970827 4.614602e+06 8547365.9
## 303 1215.997308 4.645314e+06 8560161.4
## 304 1220.023789 4.676129e+06 8560161.4
## 305 1224.050270 4.707045e+06 8624138.7
## 306 1228.076752 4.738064e+06 8636934.1
## 307 1232.103233 4.769184e+06 8662525.1
## 308 1236.129714 4.800406e+06 8916301.6
## 309 1240.156195 4.831730e+06 8929452.5
## 310 1244.182676 4.863156e+06 8929452.5
## 311 1248.209157 4.894683e+06 8929452.5
## 312 1252.235638 4.926313e+06 8955754.3
## 313 1256.262120 4.958044e+06 8982056.1
## 314 1260.288601 4.989878e+06 9087263.2
## 315 1264.315082 5.021813e+06 9087263.2
## 316 1268.341563 5.053850e+06 9113565.0
## 317 1272.368044 5.085989e+06 9192470.3
## 318 1276.394525 5.118229e+06 9205621.2
## 319 1280.421006 5.150572e+06 9258224.7
## 320 1284.447488 5.183016e+06 9258224.7
## 321 1288.473969 5.215563e+06 9271375.6
## 322 1292.500450 5.248211e+06 9271375.6
## 323 1296.526931 5.280961e+06 9337130.0
## 324 1300.553412 5.313813e+06 9416035.4
## 325 1304.579893 5.346767e+06 9429186.2
## 326 1308.606375 5.379822e+06 9441585.7
## 327 1312.632856 5.412980e+06 9441585.7
## 328 1316.659337 5.446239e+06 9495692.2
## 329 1320.685818 5.479601e+06 9566110.2
## 330 1324.712299 5.513064e+06 9566110.2
## 331 1328.738780 5.546629e+06 9566110.2
## 332 1332.765261 5.580296e+06 9566110.2
## 333 1336.791743 5.614064e+06 9635732.5
## 334 1340.818224 5.647935e+06 9663581.5
## 335 1344.844705 5.681907e+06 9663581.5
## 336 1348.871186 5.715982e+06 9669488.8
## 337 1352.897667 5.750158e+06 9669488.8
## 338 1356.924148 5.784436e+06 9698181.7
## 339 1360.950630 5.818816e+06 9726874.5
## 340 1364.977111 5.853298e+06 9741220.9
## 341 1369.003592 5.887881e+06 9741220.9
## 342 1373.030073 5.922567e+06 9784260.2
## 343 1377.056554 5.957354e+06 9812953.1
## 344 1381.083035 5.992244e+06 9850438.9
## 345 1385.109516 6.027235e+06 9850438.9
## 346 1389.135998 6.062328e+06 9942070.9
## 347 1393.162479 6.097522e+06 9987886.9
## 348 1397.188960 6.132819e+06 10003158.8
## 349 1401.215441 6.168218e+06 10048974.8
## 350 1405.241922 6.203718e+06 10064246.8
## 351 1409.268403 6.239321e+06 10099881.5
## 352 1413.294885 6.275025e+06 10099881.5
## 353 1417.321366 6.310831e+06 10154299.0
## 354 1421.347847 6.346739e+06 10195695.1
## 355 1425.374328 6.382749e+06 10195695.1
## 356 1429.400809 6.418860e+06 10246420.0
## 357 1433.427290 6.455074e+06 10263328.2
## 358 1437.453771 6.491389e+06 10280236.5
## 359 1441.480253 6.527806e+06 10280236.5
## 360 1445.506734 6.564326e+06 10297144.8
## 361 1449.533215 6.600947e+06 10297144.8
## 362 1453.559696 6.637669e+06 10288040.4
## 363 1457.586177 6.674494e+06 10342667.1
## 364 1461.612658 6.711421e+06 10342667.1
## 365 1465.639140 6.748449e+06 10342667.1
## 366 1469.665621 6.785579e+06 10342667.1
## 367 1473.692102 6.822812e+06 10342667.1
## 368 1477.718583 6.860146e+06 10360876.0
## 369 1481.745064 6.897582e+06 10397293.9
## 370 1485.771545 6.935119e+06 10415502.8
## 371 1489.798026 6.972759e+06 10451920.6
## 372 1493.824508 7.010501e+06 10470129.6
## 373 1497.850989 7.048344e+06 10621342.8
## 374 1501.877470 7.086289e+06 10621342.8
## 375 1505.903951 7.124336e+06 10621342.8
## 376 1509.930432 7.162485e+06 10683094.8
## 377 1513.956913 7.200736e+06 10703678.8
## 378 1517.983394 7.239089e+06 10703678.8
## 379 1522.009876 7.277544e+06 10703678.8
## 380 1526.036357 7.316100e+06 10759816.9
## 381 1530.062838 7.354758e+06 10759816.9
## 382 1534.089319 7.393519e+06 10753668.5
## 383 1538.115800 7.432381e+06 10821301.6
## 384 1542.142281 7.471344e+06 10843846.0
## 385 1546.168763 7.510410e+06 10843846.0
## 386 1550.195244 7.549578e+06 10843846.0
## 387 1554.221725 7.588847e+06 10843846.0
## 388 1558.248206 7.628219e+06 10934023.5
## 389 1562.274687 7.667692e+06 10934023.5
## 390 1566.301168 7.707267e+06 10934023.5
## 391 1570.327649 7.746944e+06 10934023.5
## 392 1574.354131 7.786723e+06 11362366.7
## 393 1578.380612 7.826604e+06 11409709.9
## 394 1582.407093 7.866586e+06 11486954.0
## 395 1586.433574 7.906671e+06 11486954.0
## 396 1590.460055 7.946857e+06 11599082.7
## 397 1594.486536 7.987145e+06 11651686.2
## 398 1598.513018 8.027535e+06 11677988.0
## 399 1602.539499 8.068027e+06 11677988.0
## 400 1606.565980 8.108621e+06 11677988.0
## 401 1610.592461 8.149317e+06 11677988.0
## 402 1614.618942 8.190114e+06 11704289.8
## 403 1618.645423 8.231014e+06 11704289.8
## 404 1622.671904 8.272015e+06 11730591.5
## 405 1626.698386 8.313118e+06 11756893.3
## 406 1630.724867 8.354323e+06 11835798.6
## 407 1634.751348 8.395630e+06 11835798.6
## 408 1638.777829 8.437039e+06 11914704.0
## 409 1642.804310 8.478549e+06 11914704.0
## 410 1646.830791 8.520162e+06 11914704.0
## 411 1650.857273 8.561876e+06 11941005.7
## 412 1654.883754 8.603692e+06 11967307.5
## 413 1658.910235 8.645610e+06 11967307.5
## 414 1662.936716 8.687630e+06 11967307.5
## 415 1666.963197 8.729752e+06 11993609.3
## 416 1670.989678 8.771975e+06 11993609.3
## 417 1675.016159 8.814301e+06 11993609.3
## 418 1679.042641 8.856728e+06 11993609.3
## 419 1683.069122 8.899258e+06 11993609.3
## 420 1687.095603 8.941889e+06 11993609.3
## 421 1691.122084 8.984622e+06 12019911.1
## 422 1695.148565 9.027457e+06 12072514.6
## 423 1699.175046 9.070393e+06 12072514.6
## 424 1703.201528 9.113432e+06 12072514.6
## 425 1707.228009 9.156572e+06 12072514.6
## 426 1711.254490 9.199815e+06 12098816.4
## 427 1715.280971 9.243159e+06 12098816.4
## 428 1719.307452 9.286605e+06 12125118.2
## 429 1723.333933 9.330153e+06 12125118.2
## 430 1727.360414 9.373802e+06 12125118.2
## 431 1731.386896 9.417554e+06 12125118.2
## 432 1735.413377 9.461408e+06 12125118.2
## 433 1739.439858 9.505363e+06 12125118.2
## 434 1743.466339 9.549420e+06 12125118.2
## 435 1747.492820 9.593579e+06 12125118.2
## 436 1751.519301 9.637840e+06 12125118.2
## 437 1755.545782 9.682203e+06 12125118.2
## 438 1759.572264 9.726668e+06 12125118.2
## 439 1763.598745 9.771234e+06 12125118.2
## 440 1767.625226 9.815903e+06 12151419.9
## 441 1771.651707 9.860673e+06 12151419.9
## 442 1775.678188 9.905545e+06 12177721.7
## 443 1779.704669 9.950519e+06 12177721.7
## 444 1783.731151 9.995595e+06 12058590.1
## 445 1787.757632 1.004077e+07 12058590.1
## 446 1791.784113 1.008605e+07 12058590.1
## 447 1795.810594 1.013143e+07 12058590.1
## 448 1799.837075 1.017692e+07 12058590.1
## 449 1803.863556 1.022250e+07 12058590.1
## 450 1807.890037 1.026819e+07 12058590.1
## 451 1811.916519 1.031398e+07 12058590.1
## 452 1815.943000 1.035987e+07 12058590.1
## 453 1819.969481 1.040586e+07 12058590.1
## 454 1823.995962 1.045196e+07 12114288.0
## 455 1828.022443 1.049815e+07 12197834.8
## 456 1832.048924 1.054445e+07 12197834.8
## 457 1836.075406 1.059085e+07 12197834.8
## 458 1840.101887 1.063735e+07 12197834.8
## 459 1844.128368 1.068396e+07 12197834.8
## 460 1848.154849 1.073066e+07 12225683.8
## 461 1852.181330 1.077747e+07 12225683.8
## 462 1856.207811 1.082438e+07 12225683.8
## 463 1860.234292 1.087139e+07 12225683.8
## 464 1864.260774 1.091851e+07 12225683.8
## 465 1868.287255 1.096572e+07 12250051.6
## 466 1872.313736 1.101304e+07 12250051.6
## 467 1876.340217 1.106046e+07 12279641.1
## 468 1880.366698 1.110798e+07 12279641.1
## 469 1884.393179 1.115560e+07 12279641.1
## 470 1888.419661 1.120332e+07 12338820.1
## 471 1892.446142 1.125115e+07 12467041.2
## 472 1896.472623 1.129908e+07 12530165.5
## 473 1900.499104 1.134711e+07 12530165.5
## 474 1904.525585 1.139524e+07 12530165.5
## 475 1908.552066 1.144347e+07 12512130.0
## 476 1912.578547 1.149181e+07 12545946.6
## 477 1916.605029 1.154025e+07 12579763.1
## 478 1920.631510 1.158879e+07 12613579.7
## 479 1924.657991 1.163743e+07 12613579.7
## 480 1928.684472 1.168617e+07 12681212.8
## 481 1932.710953 1.173502e+07 12681212.8
## 482 1936.737434 1.178396e+07 12681212.8
## 483 1940.763916 1.183301e+07 12715029.4
## 484 1944.790397 1.188216e+07 12673409.0
## 485 1948.816878 1.193141e+07 12709826.8
## 486 1952.843359 1.198077e+07 12709826.8
## 487 1956.869840 1.203023e+07 12709826.8
## 488 1960.896321 1.207978e+07 12782662.5
## 489 1964.922802 1.212944e+07 12610505.5
## 490 1968.949284 1.217920e+07 12610505.5
## 491 1972.975765 1.222907e+07 12610505.5
## 492 1977.002246 1.227903e+07 12545946.6
## 493 1981.028727 1.232910e+07 12545946.6
## 494 1985.055208 1.237927e+07 12545946.6
## 495 1989.081689 1.242954e+07 12593289.8
## 496 1993.108170 1.247991e+07 12593289.8
## 497 1997.134652 1.253039e+07 12593289.8
## 498 2001.161133 1.258097e+07 12593289.8
## 499 2005.187614 1.263164e+07 12677455.4
## 500 2009.214095 1.268243e+07 12677455.4
## 501 2013.240576 1.273331e+07 12677455.4
## 502 2017.267057 1.278429e+07 12677455.4
## 503 2021.293539 1.283538e+07 12677455.4
## 504 2025.320020 1.288657e+07 12677455.4
## 505 2029.346501 1.293786e+07 12677455.4
## 506 2033.372982 1.298925e+07 12677455.4
## 507 2037.399463 1.304074e+07 12677455.4
## 508 2041.425944 1.309234e+07 12677455.4
## 509 2045.452425 1.314403e+07 12730059.0
## 510 2049.478907 1.319583e+07 12730059.0
## 511 2053.505388 1.324773e+07 12730059.0
## 512 2057.531869 1.329974e+07 12730059.0
## 513 2061.558350 1.335184e+07 12730059.0
There are a lot of elements i the plot of K that can be explained. First, the Kpois(r) line in Red is the theoretical value of K for each distance window (r) under a Poisson assumption of Complete Spatial Randomness.[Practical book] The Black line is the estimated values of K accounting for the effects of the edge of the study area.[Practical book]
(Here, the correction specifies how points towards the edge are dealt with, in this case, border means that points towards the edge are ignored for the calculation but are included for the central points.)
If the value of K falls above the line, it means that there is cluster of the data at the distance. If the value of K is below the line, there shows a dispersion of the data. From the graph, we can see that up until distances of around 2000 metres, the location of the eviction seems to be clustered in Harrow. (however, at around 1500 m, the distribution appears random and then dispersed between about 1600 and 2100 metres.)
# Load the package
library(raster)
##
## Attaching package: 'raster'
## The following object is masked from 'package:janitor':
##
## crosstab
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
## The following objects are masked from 'package:MASS':
##
## area, select
## The following object is masked from 'package:nlme':
##
## getData
## The following objects are masked from 'package:spatstat.geom':
##
## area, rotate, shift
library(fpc)
#first check the coordinate reference system of the CD502 spatial polygon:
st_geometry(NYDistricts)
## Geometry set for 71 features
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 563069.7 ymin: 4483098 xmax: 609762.3 ymax: 4529952
## Projected CRS: WGS 84 / UTM zone 18N
## First 5 geometries:
## MULTIPOLYGON (((595104.1 4522026, 595098.5 4521...
## MULTIPOLYGON (((597308.6 4510424, 597266.4 4510...
## MULTIPOLYGON (((594364.8 4521323, 594363 452132...
## MULTIPOLYGON (((593253.9 4504080, 593249.3 4504...
## MULTIPOLYGON (((593432.1 4523996, 593429 452399...
DBSCAN requires to input two parameters: 1. Epsilon - this is the radius within which the algorithm with search for clusters 2. MinPts - this is the minimum number of points that should be considered a cluster
Based on the results of the Ripley’s K analysis earlier, we can see that we are getting clustering up to a radius of around 2000m, with the largest bulge in the graph at around 1700m. Therefore, 1700m is probably a good place to start and we will begin by searching for clusters of at least 4 points…
#first extract the points from the spatial points data frame
NewEvictionsCDPoints <- NewEvictionsCD %>%
coordinates(.)%>%
as.data.frame()
#now run the dbscan analysis
db <- NewEvictionsCDPoints %>%
fpc::dbscan(.,eps = 800, MinPts = 6)
#now plot the result
plot(db, NewEvictionsCDPoints, main = "DBSCAN Output", frame = F)
plot(NYDistricts$geometry, add=T)
Use kNNdistplot() to find a suitable eps value based on the ‘knee’ in the plot…
# k is no of nearest neighbours used, use min points
library(dbscan)
##
## Attaching package: 'dbscan'
## The following object is masked from 'package:fpc':
##
## dbscan
NewEvictionsCDPoints%>%
dbscan::kNNdistplot(.,k=6)
https://www.datanovia.com/en/lessons/dbscan-density-based-clustering-essentials/
Advanced DBSCAN
library(ggplot2)
db
## dbscan Pts=116 MinPts=6 eps=800
## 0 1 2 3 4
## border 22 3 1 4 0
## seed 0 21 25 32 8
## total 22 24 26 36 8
db$cluster
## [1] 0 3 0 1 1 0 2 3 0 3 2 2 3 2 3 4 3 3 3 3 1 3 3 3 2 2 3 3 0 3 1 1 4 1 0 3 3
## [38] 0 2 2 2 0 2 3 2 0 1 0 3 3 3 1 3 2 3 3 2 3 3 0 2 1 1 1 0 3 2 2 2 4 3 1 1 0
## [75] 2 4 0 1 2 0 0 3 1 2 3 1 4 3 4 2 1 1 3 1 0 3 2 0 1 2 2 1 0 4 1 0 0 3 3 4 2
## [112] 3 1 2 0 1
# add this cluster membership info back into our dataframe
NewEvictionsCDPoints <- NewEvictionsCDPoints %>%
mutate(dbcluster=db$cluster)
# Next we are going to create some convex hull polygons to wrap around the points in our clusters.
chulls <- NewEvictionsCDPoints %>%
group_by(dbcluster) %>%
dplyr::mutate(hull = 1:n(),
hull = factor(hull, chull(coords.x1, coords.x2)))%>%
arrange(hull)
#chulls2 <- ddply(BluePlaquesSubPoints, .(dbcluster),
# function(df) df[chull(df$coords.x1, df$coords.x2), ])
# As 0 isn’t actually a cluster (it’s all points that aren’t in a cluster) drop it from the dataframe
chulls <- chulls %>%
filter(dbcluster >=1)
# ggplot
dbplot <- ggplot(data = NewEvictionsCDPoints,
aes(coords.x1,coords.x2, colour=dbcluster, fill=dbcluster))
#add the points in
dbplot <- dbplot + geom_point()
#now the convex hulls
dbplot <- dbplot + geom_polygon(data = chulls,
aes(coords.x1,coords.x2, group=dbcluster),
alpha = 0.5)
#now plot, setting the coordinates to scale correctly and as a black and white plot
#(just for the hell of it)...
dbplot + theme_bw() + coord_equal()
Nicer way of DBSCAN
#convex hulls to wrap around points
chulls2 <- data.frame()
for (cluster in 1:max(NewEvictionsCDPoints$dbcluster)) {
cluster_data <- NewEvictionsCDPoints %>%
filter(dbcluster == cluster)
ch <- chull(cluster_data$coords.x1, cluster_data$coords.x2)
chulls2 <- chulls2 %>%
bind_rows(cluster_data[c(ch), ])
}
# ggplot
dbplot <- ggplot(data = NewEvictionsCDPoints,
aes(coords.x1,coords.x2, colour=dbcluster, fill=dbcluster))
#add the points in
dbplot <- dbplot + geom_point()
#now the convex hulls
dbplot <- dbplot + geom_polygon(data = chulls2,
aes(coords.x1,coords.x2, group=dbcluster),
alpha = 0.5)
#now plot, setting the coordinates to scale correctly and as a black and white plot
#(just for the hell of it)...
dbplot + theme_bw() + coord_equal()
library(here)
library(janitor)
library(dplyr)
# join the New York District data and eviction data
#Joinfun <- function(data1, data2){
#output<- data1%>%
# st_join(NYDistricts,.)%>%
# add_count(boro_cd, name="eviction-borough")
# return(output)
#}
#NYEvictions <- Joinfun(NewEvictions, NYDistricts)
#NYEvictions
#summary(NYEvictions)
#NewEvictions <- NewEvictions[NYEvictions,]
#tm_shape(NYDistricts) +
# tm_polygons(col = NA, alpha = 0.5) +
#tm_shape(NYEvictions) +
# tm_dots(col = "blue")
library(sf)
points_sf_joined <- NYDistricts %>%
st_join(NewEvictions) %>%
add_count(boro_cd) %>%
janitor::clean_names() %>%
#calculate area
mutate(area=st_area(.)) %>%
#then density of the points per ward
mutate(density=n/area) %>%
#select density and some other variables
dplyr::select(density, boro_cd, n)
points_sf_joined <- points_sf_joined %>%
group_by(boro_cd) %>%
summarise(density = first(density),
boro_cd = first(boro_cd),
plaquecount= first(n))
tm_shape(points_sf_joined) +
tm_polygons("density",
style="jenks",
palette="PuOr",
midpoint=NA,
popup.vars=c("boro_cd", "density"),
title="Evictions Density")
So, from the map, it looks as though we might have some clustering of blue plaques in the Manhatten so let’s check this with Moran’s I and some other statistics.
Before being able to calculate Moran’s I and any similar statistics, we need to first define a W_ij spatial weights matrix.
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
#First calculate the centroids of all Wards in New York City
coordsW <- points_sf_joined%>%
st_centroid()%>%
st_geometry()
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
plot(coordsW,axes=TRUE)
#create a neighbours list
LWard_nb <- points_sf_joined %>%
poly2nb(., queen=T)
summary(LWard_nb)
## Neighbour list object:
## Number of regions: 71
## Number of nonzero links: 314
## Percentage nonzero weights: 6.228923
## Average number of links: 4.422535
## Link number distribution:
##
## 1 2 3 4 5 6 7 8
## 1 8 11 16 18 10 6 1
## 1 least connected region:
## 62 with 1 link
## 1 most connected region:
## 33 with 8 links
Average number of links is 4.422535
#plot the neighbours
plot(LWard_nb, st_geometry(coordsW), col="red")
#add a map underneath
plot(points_sf_joined$geometry, add=T)
#create a spatial weights matrix from these weights
#Lward.lw <- LWard_nb %>%
# nb2mat(., style="B")
#sum(Lward.lw)
#summary(Lward.lw)
Now we have defined our W_ij matrix, we can calculate the Moran’s I and other associated statistics. However, Moran’s I requires a spatial weight list type object as opposed to matrix.
Lward.lw <- LWard_nb %>%
nb2listw(., style="C")
Moran’s I test tells us whether we have clustered values (close to 1) or dispersed values (close to -1)
I_LWard_Global_Density <- points_sf_joined %>%
pull(density) %>%
as.vector()%>%
moran.test(., Lward.lw)
I_LWard_Global_Density
##
## Moran I test under randomisation
##
## data: .
## weights: Lward.lw
##
## Moran I statistic standard deviate = 8.238, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.593759367 -0.014285714 0.005447915
This tells us whether similar values or dissimilar values are clustering
C_LWard_Global_Density <-
points_sf_joined %>%
pull(density) %>%
as.vector()%>%
geary.test(., Lward.lw)
C_LWard_Global_Density
##
## Geary C test under randomisation
##
## data: .
## weights: Lward.lw
##
## Geary C statistic standard deviate = 4.0732, p-value = 2.319e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.46436106 1.00000000 0.01729311
This tells us whether high or low values are clustering. If G > Expected = High values clustering; if G < expected = low values clustering
G_LWard_Global_Density <-
points_sf_joined %>%
pull(density) %>%
as.vector()%>%
globalG.test(., Lward.lw)
G_LWard_Global_Density
##
## Getis-Ord global G statistic
##
## data: .
## weights: Lward.lw
##
## standard deviate = 7.2857, p-value = 1.6e-13
## alternative hypothesis: greater
## sample estimates:
## Global G statistic Expectation Variance
## 3.048722e-02 1.428571e-02 4.945053e-06
So the global statistics are indicating that we have spatial autocorrelation of Blue Plaques in London:
The Moran’s I statistic = 0.59 (remember 1 = clustered, 0 = no pattern, -1 = dispersed) which shows that we have some distinctive clustering
The Geary’s C statistic = 0.46 (remember Geary’s C falls between 0 and 2; 1 means no spatial autocorrelation, <1 - positive spatial autocorrelation or similar values clustering, >1 - negative spatial autocorreation or dissimilar values clustering) which shows that similar values are clustering
The General G statistic = G > expected, so high values are tending to cluster.
We can now also calculate local versions of the Moran’s I statistic (for each Ward) and a Getis Ord G*_i statistic to see where we have hot-spots…
#use the localmoran function to generate I for each ward in the city
I_LWard_Local_count <- points_sf_joined %>%
pull(plaquecount) %>%
as.vector()%>%
localmoran(., Lward.lw)%>%
as_tibble()
I_LWard_Local_Density <- points_sf_joined %>%
pull(density) %>%
as.vector()%>%
localmoran(., Lward.lw)%>%
as_tibble()
#what does the output (the localMoran object) look like?
slice_head(I_LWard_Local_Density, n=5)
## # A tibble: 5 × 5
## Ii E.Ii Var.Ii Z.Ii `Pr(z != E(Ii))`
## <localmrn> <localmrn> <localmrn> <localmrn> <localmrn>
## 1 -0.001116426 -0.0015044943 0.02372412 0.002519496 0.9979897
## 2 -0.007248032 -0.0007965765 0.01222156 -0.058357190 0.9534641
## 3 -0.056804476 -0.0008635776 0.01324855 -0.486010031 0.6269601
## 4 -0.158916086 -0.0041275381 0.06309378 -0.616233868 0.5377402
## 5 0.011812369 -0.0008714819 0.01296872 0.111378795 0.9113160
points_sf_joined <- points_sf_joined %>%
mutate(plaque_count_I = as.numeric(I_LWard_Local_count$Ii))%>%
mutate(plaque_count_Iz =as.numeric(I_LWard_Local_count$Z.Ii))%>%
mutate(density_I =as.numeric(I_LWard_Local_Density$Ii))%>%
mutate(density_Iz =as.numeric(I_LWard_Local_Density$Z.Ii))
We’ll set the breaks manually based on the rule that data points >2.58 or <-2.58 standard deviations away from the mean are significant at the 99% level (<1% chance that autocorrelation not present); >1.96 - <2.58 or <-1.96 to >-2.58 standard deviations are significant at the 95% level (<5% change that autocorrelation not present).
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
MoranColours<- rev(brewer.pal(8, "RdGy"))
tm_shape(points_sf_joined) +
tm_polygons("plaque_count_Iz",
style="fixed",
breaks=breaks1,
palette=MoranColours,
midpoint=NA,
title="Local Moran's I, Blue Plaques in London")
This map shows some areas in the centre of London that have relatively high scores, indicating areas with lots of blue plaques neighbouring other areas with lots of blue plaques.
the Getis Ord G*_i statisic for hot and cold spots
Gi_LWard_Local_Density <- points_sf_joined %>%
pull(density) %>%
as.vector()%>%
localG(., Lward.lw)
head(Gi_LWard_Local_Density)
## [1] -0.002519496 0.058357190 -0.486010031 -0.616233868 -0.111378795
## [6] -0.479051443
points_sf_joined <- points_sf_joined %>%
mutate(density_G = as.numeric(Gi_LWard_Local_Density))
GIColours<- rev(brewer.pal(8, "RdBu"))
#now plot on an interactive map
tm_shape(points_sf_joined) +
tm_polygons("density_G",
style="fixed",
breaks=breaks1,
palette=GIColours,
midpoint=NA,
title="Gi*, Blue Plaques in London")
#### Other variables can be checked in the csv file.